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O Abstract : 
O 

We investigate a reformulation of the dynamics of interacting fermion systems in 
<^ terms of a stochastic extension of Time Dependent Hartree-Fock equations. The noise is 
CN found from a path-integral representation of the evolution operator and allows to interpret the 
exact N-body state as a coherent average over Slater determinants evolving under the random 

O 

^ mean-fied. The full density operator and the expectation value of any observable are then 
O reconstructed using pairs of stochastic uncorrelated wave functions. The imaginary time 
propagation is also presented and gives a similar stochastic one-body scheme which 
converges to the exact ground state without developing a sign problem. In addition, the 
. ^ growth of statistical errors is examined to show that the stochastic formulation never explode 
in a finite dimensional one-body space. Finally, we consider initially correlated systems and 
present some numerical implementations in exactly soluble models to analyse the precision 
and the stability of the approach in practical cases. 
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I - Introduction 



Even with present-day computing facilities, the theoretical study of the structure and 
the dynamics of many-fermion systems, like nuclei, atomic clusters or quantum dots, remains 
a formidable task which requires approximations to step down to feasibility. In such a context, 
Hartree-Fock (HF) theory [1] is usually considered as the basic tool. The complex many -body 
problem is then reduced to an effective single-particle description in which the interaction is 
communicated through a common and self consistent mean-field. In nuclear physics, static 
HF calculations provide a very good starting point to investigate many properties of the 
ground state [2]. The time-dependent model (TDHF) is also successful at low energy in 
heavy-ion reactions [2] and for the nonlinear electron dynamics in metal clusters [3]. 

Nevertheless, the HF description only treats in average the two-body potential and 
thus needs to be improved for strongly correlated systems where the residual interaction can 
generate coherent superpositions of many mean-fields and gives rise to collision-type 
processes in the dynamical evolution. In consequence, various schemes have been developed 
to take into account correlation effects beyond the HF approximation. For example, in the 
condensed matter problem and in nuclear spectroscopy, recent developments make use of a 
path integral representation of the Boltzmann operator to select, in a stochastic way, mean- 
fields dominating the structure of low-lying states [4,5]. Concerning dynamical problems, 
extensions of TDHF are also possible by adding a Boltzmann-like collision term deduced 
from an appropriate truncation in the hierarchy of equations describing the time-evolution of 
reduced density-matrix [6]. However, realistic applications of such a quantum transport 
theory remain difficult due to the large number of rarely occupied orbitals needed in the 
simulation [7]. Most of collisional mean-field calculations have thus been carried out in the 
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semi-classical limit via an integro-differential equation for the one-body phase space 
distribution, known as the Boltzmann (or Vlasov)-Uhling-Uhlenbeck (BUU or VUU) model 
[8]. In nuclear physics as in the context of metal clusters, many interesting studies of medium- 
energy processes have been done with such a transport equation concerning, for example, 
thermalization induced by binary collisions [9] or energy tranfers [10]. Nevertheless, the 
introduction of the Pauli blocked collision term becomes insufficient when fluctuations are 
crucial such as in the formation of intermediate mass fragments in heavy-ion reactions [11]. 
To investigate these situations, the BUU scheme can be improved with the help of a stochastic 
force, leading to the Boltzmann-Langevin equation (BLE) [12] which has given very 
promising results in the study of a possible connection between multifragmentation 
phenomena and the spinodal decomposition of nuclear matter [13]. In addition, the BLE 
model can be derived in the context of various quantal frameworks [14] and in particular from 
a stochastic time-dependent Hartree Fock theory [15]. After a sufficiently large time interval, 
this approach assumes that the full density can be reduced to a statistical mixing of Slater 
determinants moving in their own mean-field, each of them having a probability given by the 
residual interaction via Fermi's Golden rule. Such a scheme has only been applied in a BUU- 
type evaluation of the jumps between TDHF trajectories [16] and its validity remains to be 
shown in practical cases. 

The aim of this Letter is to present a new formulation of the many-fermion problem in 
terms of mean-field equations with a 1 particle - 1 hole (lp-lh) gaussian white noise. The 
stochastic one-body process thus constructed will have an average identical to the exact 
solution and can be used for spectroscopic or dynamical studies as illustrated in exactly 
soluble models. Statistical errors are moreover analytically calculated to show the stability of 
the proposed scheme. 
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II- From the many-body problem to an exact stochastic mean-field formulation 

To begin, let us consider the building of correlations on a Slater state ]*¥) during a 
small time-step At and under a general Hamiltonian H = ^T tj a[a j - — EK#/ a i r a + j a k a l with 

ij ^ ijkl 

a two-body interaction. If \a) denotes in general the hole states of \ X H) and \a) the particle 
orbitals, it follows immediately from the Schrodinger equation that : 



(. . At „rn^ 



At- 



a,a 2 a,a 2 



where t/(Ar) is the evolution operator and |a) the biorthogonal hole basis which satisfies 
(ctj |a 2 ) = 5„ and that have been introduced to take into account a possible non- 
orthogonality of the occupied states in the considered Slater determinant 1^). In addition, the 
coefficients in the expansion (1) are given by : 

" <0) =1 E(«|i> T, 01«> + ^E E(«il f ) («2 W ^ ««.) (% 2 ) (2-a) 

a ij ^ a ] a 2 ijkl 



Furthermore, given any two-body interaction, it is always possible to find one body operators 
Q s so that the 2p-2h amplitude (2.c) is written as : 

41^ = ^E h ^ («. I al«i> («2 1 e> 2 > (3) 

For example, a possible set {<2 s ,ft) s } can be obtained from the eigenvalue equation 
E^WW = h^^'l Q s \j) associated to the diagonalization of the hermitian matrix 

k,l 
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M {i j)( k ,i) = V ukj • Introducing X s = (Q s + Q + S )jl and P s = (Q s - Q + S )jli, the evolution (1) of the 
Slater state \¥) to first order in the time- step At becomes : 



C/(A0 |*) = Expl ^ ) Exp ^£ JXL, < a &i 

V zn / y zn a, a, 



At- 



\~[Exp 



icoAt 



Exp 



(4) 



The dynamics can then be linearized with the help of a Hubbard-Stratonovitch transformation 
[17] allowing to interpret each evolution under the square of a lp-lh operator as an infinite 
superposition of one-body evolutions, each of them in a fluctuating auxilliary field o 
distributed on a Gaussian weight. This transformation is in fact currently used in Monte-Carlo 
calculations of the equilibrium state of interacting-fermions systems at finite temperature [4]. 
However its direct application to real time problems leads to a complex distribution of the 
variables o [18] and we suggest to use as an alternative the following identity : 



Exp{i x0 2 )= ^1— J do Exp ,&cp[(l + i sign(x)) *J\x\ o oj 



f o 2 ^ 



(5) 



where x is real and O an operator. Defining the vector 6 of all the fields [o x , o P ) that are 
introduced by the path-integral (5) to linearize the evolution (4), and introducing the 

, we finally obtain : 



j~ j~ d(j d(j 

shorthand notation do G(o) = Y\ — ~ _ ~ Exp 



In 



v 



J 



U(At) \V) = Expl ^-H^ ) \doG{o) Exp £ £ AZ^ (5) a+ 



V ih 



(6) 



where 



.„ At TT(l) l + isign(co) j—, — r ,_ , , 

j 0,=X s ,P s ^ 



(7) 
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But, according to Thouless's theorem [1], the evolution of the Slater determinant |*P) under 
any of the fluctuating lp-lh Hamiltonians of the auxiliary field representation (6) gives a 
neighboring pure independent-particle vector whose hole states are linear combinations of the 
initial particle and hole orbitals. To first-order in At, the correlated wave function (6) can thus 
be brought, with the help of (2. a) and (2.b), in the form : 

U(At)\V) = \daG(a)Ylal Aa{ r } \ ) (8) 

a 

with 

|Aa,(<T)) = E |«2> AF^ fc)AZ^(<j), (9) 
the deterministic contribution AF being defined by : 

L(«i|i) T ij (i| a 2 ) + ^L E(«i|0 ( & 3 \j) v mi ( k \ a 2) { l \ a -: 



ih 



a l a 2 



2' 

a 3 ijkl 



(10) 



Together with the expression (7) of the <7 -dependent element AZ, and after some algebraic 
manipulations, the variation (9) of the holes states can in fact also be written as : 

1 + i sign(a> s ) 



\Aa(a)) = - 

1 v ih 



T+V(p)-^pV(p) 



\a 



A^ + E £ %-^o 0i $^t{l-p)0,\a) 



s O s =X s ,P, 



(ID 

where p = ^\a) (a\ is the one-body density and V(p) the mean-field potential with matrix 

a 

elements (i| V(p) \j) = ^V ikjl p lk . In the end, the full dynamics of an uncorrelated state 

u 

during a small time-step At under a one and two-body Hamiltonian, can be represented as 
the coherent average of Slater determinants that have evolved with a mean-field Hamiltonian 
supplemented with a fluctuating lp-lh potential. This noise, which comes from the 2p-2h 
residual interaction, is in addition non-hermitian and depends linearly on external number- 
fields, each distributed on a unit gaussian weight. We also emphasize that the deterministic 

part differs from the standard Hartree-Fock approach by including the term - ^pV(p) which 



-6- 



in fact arises naturally to take into account the difference beween the total mean-field energy 
and the sum of the HF eigenvalues [2]. Last but by no means least, the formulation (8)-(l 1) of 
the dynamics can be extended to many time-steps leading to a wave function at any time 
t = NAt that is always the expectation of random Slater determinants. The Fermi sea 
associated to these states is then evaluated by a recursive application of (11) with a different 
stochastic vector field 6 at each iteration. Such an algorithm corresponds in fact to the 
discretized realization of a Markovian Langevin equation in the one-body Hilbert space in 
accordance to the strong Euler scheme [19]. In consequence, going to the continuous limit 
At — > , and denoting by E (L ) the expectation of a random functional, we can state that : 

f \ 



u( t )\v)=E n<, 



) 

V a 



(12) 



where the evolution of the hole states \cc(t)) is determined by the following Ito stochastic 
differential equation [19] with the initial condition |a(0)) = \oc): 

d\a(t)) = ± T + v{p(t))-^p(t)V{p(t)) \a(t))dt 

+ 1 I 1 + - T K) VKf(l- Pit)) O, \a(t)) dW 0< (t) (13) 

where W refer to independent real Wiener processes with vanishing enesemble averages and 
that obey to the Ito stochastic calculus [20] : 

E(dW (t)) = Q (14.a) 
dW (t)dW (t) = dt 8 (14.b) 

It is also important to specify that such a formulation provides a reinterpretation of the exact 
evolution of the N-body density operator D{t) in terms of a mean over dyadics in which the 
bra and the ket are different Slater determinants fj<z„( ( )| ), \\ a+ a:(t)\ ) w ^ m an d | Of' (/■)) 

a a' 

evolving independently under the stochastic Hartree-Fock equation (12) : 
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D(t) = U(t) |V) (V\ U + {t) = E ( a + aAt) \ ) ( | a av(0 ... a a , M ) (15) 

Similarly, the expectation value of any observable A will be expressed by : 

(A)(t) = Tr[D(t)A] = E(( \ a^ (t) ... a a , M A | ) ) (16) 

The necessary use of pairs of stochastic uncorrelated wave functions to reconstruct the exact 
full density operator consitutes in fact the originality of the proposed approach in comparison 
with previous schemes [15] where a quantum incoherent propagation of states was assumed. 
By a totally different way and for interacting fermions systems, we thus end up with the same 
representation of the N-body density matrix as those obtained for the dynamics of Bose- 
Einstein condensates [21]. 



Finally, the stochastic mean-field interpretation (12)-(14) of the evolution process can 
be extended to imaginary time to obtain the ground state of the system. The only modification 
in the previous demonstration is now the use of the standard Hubbard- Stratonovitch 
transformation [17] which leads to the following representation of the Boltzmann operator 
when it acts on a Slater determinant | *F) = Y\ a a | ) : 



Exp(-f3H)\V) = E(Ha^ 

V a ) 



(17) 



with |a(0)) = | a) and 



d\a(p)) = \r + v{p(p))-^p(p) V(p((3)) 



HP)) dp 



+ I I 

s (> \ ./' 



(l + sign((O s )) + i (l - sign((D s )) fh 



^(l-p(P))O s \a(p))dW (p) 

(18) 



The correlated many-body ground state can then be obtained in the limit of large P where 
Exp(~PH) bahaves like a filter eliminating all the overlaps with excited eigenstates. The 
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stochastic one-body interpretation (17)-(18) is in fact particularly interesting since it never 
suffers from the sign problem which occurs in the Monte-Carlo calculations of expectations 
values in the canonical ensemble. 



Ill- Statistical fluctuations and stability of the stochastic approach 

We now examine the growth of statistical errors in the stochastic sheme (12)-(13) or 
(17)-(18). This can be achieved by calculating the average of the norm of the deviation 
between the exact evolution and a stochastic realization. For time-dependent problems, we are 
thus interesting in the indicator : x{t) = E ( (<p(t)\<p(t)) ) with \<p(t)) = U(t) |¥) - Y[K(t)\ ) 

a 

where the single-particle wave-functions \a(t)) evolve according to the Ito equation (13). 
Using Eq. (12) and the unitarity of the evolution operator, one then immediatly checks that : 

X (t)= e ((nt)\nt)))-i^h\^(t)}=u<(,) \ ) < 19 ) 

a 

Furthermore (^(/)|^(/)) corresponds to the determinant of the overlap matrix 
g a a (t) = (g^/) | C£ 2 (/)) associated to the Fermi sea of ^(t)) and which evolves according to :. 
dg ai a 2 (t) = {da 1 (t)\a 2 (t)) + (a 1 (t)\da 2 (t)) + (da 1 (t)\da 2 (t)) 

= I KK«iW| 0,(1-^0)0, wo) 

^ s O s =X s ,P s 

(20) 

where the stochastic Hartree-Fock equations (13) and Ito's differentiation rules [20] have 
been used. In consequence the variation of % during an infinitesimal time dt is given by : 

d X {t) = E( (V(f)| V(f)> Tr(dg(t) g-\t)) ) (21) 



leading finally to : 



d X (t)=*E 



V s O s =X s ,P s 



(22) 
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where (^ W ^(0 S ) refer to the quantal fluctuations of the one-body operators O s in the 
uncorrected stochastic wave function Such a result then implies that the proposed 

approach will never explode in a finite dimensional one-body space : in this case, a trivial 
upper bound for the variance <j\, (t) (O s ) is in fact given by the square of the largest eigenvalue 

X of the observable O s and thus dx(t) < — ^(0^ ^2 1^1 ^o> which after time 

2 s o s =x s ,p s 

integration gives : 



X(t) < Exp 



( . \ 

2 



j o=x s ,p s j 



o. 



(23) 



IV- Simulation examples 

As a first a step towards the use of the previous stochastic formulation in practical 
problems, we investigate some numerical implementations in exactly soluble models. More 
precisely, we consider a system of Q. fermions distributed among a number n of energy 
orbitals, each of which is Q. -degenerate. Single particle states are thus labelled by two 
numbers i = 0,L ,n - 1 for the level and (O = 1,L ,£l to take into account degeneracy within 
each orbital. In addition, the Hamiltonian of the system is expanded onto the bilinear 

£1 

operators =^« ! ^ fi) which generates a U(n) Lie algebra : usually, it includes the one- 

0) = 1 

n-1 

body contribution ^e ; . G„ of the individual levels with energies £, and a two-body 



i=0 



interaction written in terms of products of the generators G, y . This non trivial model is 

analytically soluble by group theoretical techniques in some cases and mimics in particular 
the shell-model picture of the atomic nucleus. For n = 2 levels, it corresponds to the standard 
Lipkin-Meshkov-Glick model [22] which is usually used to check the validity of approximate 
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many-body techniques. The three-orbital case has also been studied in particular in the semi- 
classical limit where the dynamics can exhibit a chaotic behavior [23]. 



We first focus on the two-level model with the following Hamiltonian : 

H = eJ z + VJ 2 x (24) 
where the quasi-spin operators J x = (G 10 + G 01 )/ 2, J = (G 10 - G 01 )/2i and 
J z = (G n - G 00 )/2 have been introduced. The system is initialized in an excited Hartree-Fock 

state that we choose as J I —=a^ m — ]=<^ m \ )■ This vector diagonalizes J x with the 

m=\ W2 V2 ' ) 



lowest eigenvalue m = -£1/2 in the symmetric representation of SU(2) characterized by an 
angular momentum quantum number j = £1/2. We then solve the stochastic Hartree-Fock 
equations (13) by a weak scheme [19] and use Eq. (16) to compute expectation values of 
some observables. The results are presented in Fig. 1 for £1 = 20 particles and a strength 
interaction 77 = £lV/e = 2 . After an average over 5.10 6 stochastic realizations, the expectation 
of the quasi-spin as well as those of the two and three-body operators J z , J 3 z are at all time in 
perfect agreement with the exact solution of the Schrodinger equation. Furthermore, even if 
we have observed that the associated error bars increase quasi-exponentially with time, they 
are less than 5 % at the final time thus showing the stability of the stochastic mean-field 
approach in the model. 

We have also considered the three-level extension by examining the dynamics of a 
coherent state : 

|^(z 1 ,z 2 )> = (l + |z 1 | 2 + |z 2 | 2 p /2 Mz 1 G 1 o+^ 2 o) K) = fl <^k±Jafe I ) (25) 
n 

where |*P ) = PJflo,<J ) corresponds to the Slater determinant with all particles in the ground 

(0=1 
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orbital. The Hamiltonian includes a two-body interaction having equal matrix elements for 
raising or lowering a pair of particles : 

i=0 L i,j=Q 

In addition we choose equidistant levels (£ = 0, £ x = £, e 2 =2e) and set 77 = ClV/e - 2 . For 
Q. = 10 fermions equally shared among the three levels (z t = z 2 - 1) at the initial time, the 
evolutions of the mean number of particles in each orbital are shown in Fig. 2. Once again, 
the stochastic mean-field scheme is in position to reproduce the exact dynamics but with more 
statistical fluctuations than in the previous analysis. This is in fact easily understood because 
of the larger number of Brownian processes in the stochastic Hartree-Fock equations (13). 

For spectroscopic aspects, the implementation of the propagation in imaginary time 
(18) has been realized in the two models (24) and (26), as shown in Fig. 3. In both cases, by an 
average over many simulations to reduce statistical fluctuations, the stochastic one-body 
scheme rapidly converges to the exact binding energy. 

Finally, to study the dynamics of an initial correlated state, we have investigated a 
coupling between the stochastic evolutions in real and imaginary time. In the two-level Lipkin 
model with the Hamiltonian (24), the system is supposed to be in its ground state at time t = 
and we suddenly change the sign of the two-body interaction. From a practical point of view, 

starting from the uncorrelated wave function |~J a[ m | ), the stochastic Hartree-Fock 

representation of the Boltzmann operator (17)-(18) is first applied up to a sufficiently large 
value of ft to obtain the correlated ground state. Each stochastic Slater determinant thus 
obtained is then propagated in real time with the new Hamiltonian according to Eq. (13). An 
average over a set of such one -body realizations allows after all to reconstruct the exact many- 
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body solution. This is confirmed by the results shown in Fig. 4 even if the control of statistical 
errors is more difficult than in the previous simulations. 

V- Conclusion 

We have presented a new theoretical framework to solve exactly the N-body 
Schrodinger equation for a fermionic system with binary interactions. In particular, we have 
shown that the full dynamics of an uncorrelated state can be formultated in terms of the 
coherent average of Slater determinants that have evolved with a stochastic Hartree-Fock 
equation. The noise has been derived from a path-integral representation of the evolution 
operator and comes from the linearization of the 2p-2h residual interaction by a modified 
Hubbard-Stratonovitch transformation. In addition, a similar scheme has been proposed for 
the imaginary time propagation thus leading to a stochastic one-body interpretation of the 
ground state. Finally, the coupling with the fluctuating mean-field approach in real time 
allows to reach the exact dynamics of a correlated wave function. Statistical errors have also 
been investigated to show that the proposed stochastic equations never admit a divergent 
realization in a finite dimensional one -body space. The numerical implementation in exactly 
soluble models has confirmed the formal results : in all the cases that we have studied, the 
stochastic Hartree-Fock method when averaged over many simulations give the correct result 
with a reasonable growth of the statistical spread. We hope to address a first application to a 
realistic problem in the near future. 
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Figure captions 



Fig. 1 : Time evolution of the expectation of some observables in the two-level Lipkin model 
with the Hamiltonian (24) for Q = 20 fermions and 77 = QV/£ - 2 . The initial state is the 
Slater determinant that diagonalizes J x with the eigenvalue m = - Q/2 in the space of angular 
momentum j - Q/2. The average over 5.10 6 stochastic one-body simulations is compared to 
the exact dynamics and the standard mean-field approach. Error bars in the stochastic scheme 
are not drawn because they would be roughly of the size of the points. 

Fig. 2 : Stochastic mean field dynamics of the three-level Lipkin model in comparaison to the 
exact and the Hartree-Fock evolution. The Hamiltonian takes the simple form (26) and the 
parameters are given by £ = 0, e 1 = £, £ 2 = 2e and 77 = QV/e - 2 . We have considered 

5.10 7 simulations for Q = 10 particles in the initial coherent state (25) with z x = z 2 = 1. 

Fig. 3 : Mean energy as function of the inverse temperature /3 in the stochastic Hartree-Fock 
simulation (17)-(18). In the two-level model, we consider the average over 10 7 realizations 
for a system of Q = 20 fermions with the Hamiltonian (24) and the parameters 
£ = l,V = -0.5. The initial state is the Slater determinant that diagonalizes J x with the 
eigenvalue m = - Q/2 in the space of angular momentum j - Q/2. For the three-orbital case, 
the Hamiltonian is given by Eq. (26) with £ Q = 0, £ x = 1, £ 2 =2, V = 0.2. The results 

corresponds to an average over 10 6 stochastic simultations for Q = 10 fermions initially in the 
ground orbital. 

Fig. 4 : Dynamics of an initial correlated wave function in the two-level Lipkin model for 
Q = 20 fermions with the Hamiltonian (24). For t<0, we suppose £ = 1,V = -0.1 and 
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reconstruct the ground state using the stochastic Hartree-Fock equations (18) until (3 - 5 and 

£1 

with the wave function Y\ <h,a | ) as a starting point. At time t - + , the interaction becomes 

suddenly V = 0.1 and we then apply the one-body evolution scheme (13) on each stochastic 
Slater determinant obtained with the propagation in imaginary time. An average over 10 8 of 
such simulations is compared to the exact dynamics. 
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Figure 1 
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Figure 2 
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Figure 3 
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Figure 4 




